USE ElliDef
USE NodeInfoDef
USE metric

USE Elliptic_Solvers 
IMPLICIT NONE
TYPE (NodeInfo), POINTER :: Info
INTEGER N,error
REAL (KIND=qPrec) :: dx,dy,dz,t0,t1
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:,:) :: ut,vt,wt
! first thing allocate structures for sparse matrixes
INTEGER i,j,k,nx,ny,nz
!INTEGER, PARAMETER :: solver=PARDISO_SOLVER
!INTEGER, PARAMETER :: solver=IT_SD
INTEGER, PARAMETER :: solver=IT_CG
!INTEGER, PARAMETER :: solver=IT_FOM

!INTEGER, PARAMETER :: solver=IT_GMRES
INTEGER, PARAMETER :: how_many_steps=5
INTEGER, PARAMETER :: PRECOND= 1 ! 1 use precond, 0 do not use precond
LOGICAL iterate


ALLOCATE(Info)





! populate Info structures
Info%mX=(/25,25,200/)

nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)

Info%nDim=3 ! 3D problem
! now define geometry of the problem
Info%xLower=(/-15000./2200.,-15000./2200.,0./)

Info%xUpper=(/20000./2200,15000./2200,1./)

Info%dX=(Info%xUpper(1:3)-Info%xLower(1:3))/Info%mX
Info%dmax=0.

! assume all Neumann conditions
Info%mthbc(1:6)=0


!now copy info to elli structure


ALLOCATE(Info%ElliInfo) ! allocate space for ellikit 
Info%ElliInfo%ndim=Info%ndim ! dimension of problem
Info%ElliInfo%mX=Info%mX ! size of grid
Info%ElliInfo%dX=Info%dX ! grid increments
Info%ElliInfo%mthbc=Info%mthbc ! definition of boundary conditions




   Call setprob(Info)
PRINT*,'Dimension of problem',Info%mX

CALL DEFINE_TOPO(Info) ! define the quantities used to calculate the metric

!=====================================
PRINT*,' start accelerated iterative solver'
!=====================================

CALL CPU_TIME(t0)
error=CREATE_OPERATORS_THIN(Info%ElliInfo,ACC=.true.)
CALL CPU_TIME(t1)
PRINT*,'time to create operators =',t1-t0
! Initially zero in the boundary conditions to be safe
DO i=1,2*Info%nDim
   Info%ElliInfo%bc(i)%p=zero

END DO

! check consistency of BCs
IF(.NOT.check(Info%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF
Info%ElliInfo%qe=zero
! assign local pointers to edge velocities and fill in with random values
u=>Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);
w=>Info%ElliInfo%qe(1:nx,1:ny,1:nz+1,3)
ALLOCATE(ut(1:ny,1:nz,2),vt(1:nx,1:nz,2),wt(1:nx,1:ny,2))



CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
CALL RANDOM_NUMBER(w);
! random values for edge velocities boundary conditions
!!$CALL RANDOM_NUMBER(ut)
!!$CALL RANDOM_NUMBER(vt)
!!$CALL RANDOM_NUMBER(wt)
ut=zero;vt=zero;wt=zero
        
!u=u/1.d4;v=v/1.d4;w=w/1.d4
! only need OP= o use IT wth PARDISO precon
CALL take_div(Info%ElliInfo,.false.)
PRINT*,'initial divergence is =',SQRT(SUM(Info%ElliInfo%div**2)/PRODUCT(Info%mX(1:Info%nDim)))


CALL CPU_TIME(t0)
i=POISSON_SOLVERS_THIN(Info%ElliInfo,INIT,ACC=.true.,SOLVER=solver)
CALL CPU_TIME(t1)
PRINT*,'time to set up solvers =',t1-t0
CALl CPU_TIME(t0)
CALL hodge_thin(Info%ElliInfo,-8,.false.,ut,vt,wt,SOLVER=solver)
CALL CPU_TIME(t1)

PRINT*,'time to solve =',t1-t0
i=POISSON_SOLVERS_THIN(Info%ElliInfo,FINISH,ACC=.true.,SOLVER=solver)

CALL take_div(Info%ElliInfo,.false.)
PRINT*,'residual divergence is =',SQRT(SUM(Info%ElliInfo%div**2)/PRODUCT(Info%mX(1:Info%nDim)))


   error=FREEUP_OPERATORS(Info%ElliInfo)



PRINT*,'!========================================='
PRINT*,' end accelerated solver'
PRINT*,'!========================================='

PRINT*,'!========================================='
PRINT*,' start direct solver'
PRINT*,'!========================================='

CALL CPU_TIME(t0)
error=CREATE_OPERATORS(Info%ElliInfo)
CALL CPU_TIME(t1)
PRINT*,'time to create operators =',t1-t0
! Initially zero in the boundary conditions to be safe
DO i=1,2*Info%nDim
   Info%ElliInfo%bc(i)%p=zero

END DO

! check consistency of BCs
IF(.NOT.check(Info%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF
Info%ElliInfo%qe=zero
! assign local pointers to edge velocities and fill in with random values
u=>Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);
w=>Info%ElliInfo%qe(1:nx,1:ny,1:nz+1,3)




CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
CALL RANDOM_NUMBER(w);
! random values for edge velocities boundary conditions
!!$CALL RANDOM_NUMBER(ut)
!!$CALL RANDOM_NUMBER(vt)
!!$CALL RANDOM_NUMBER(wt)
ut=zero;vt=zero;wt=zero
        
!u=u/1.d4;v=v/1.d4;w=w/1.d4
! only need OP= o use IT wth PARDISO precon
CALL take_div(Info%ElliInfo,.false.)
PRINT*,'initial divergence is =',SQRT(SUM(Info%ElliInfo%div**2)/PRODUCT(Info%mX(1:Info%nDim)))


CALL CPU_TIME(t0)
i=POISSON_SOLVERS(Info%ElliInfo,INIT, PARDISO_SOLVER)
CALL CPU_TIME(t1)
PRINT*,'time to set up solvers =',t1-t0
CALl CPU_TIME(t0)
CALL hodge_thick(Info%ElliInfo,1d-10,ut,vt,wt)
CALL CPU_TIME(t1)
PRINT*,'time to solve =',t1-t0
i=POISSON_SOLVERS(Info%ElliInfo,FINISH,PARDISO_SOLVER)

CALL take_div(Info%ElliInfo,.false.)
PRINT*,'residual divergence is =',SQRT(SUM(Info%ElliInfo%div**2)/PRODUCT(Info%mX(1:Info%nDim)))


   error=FREEUP_OPERATORS(Info%ElliInfo)

!=================================
PRINT*,'finish direct solver'
!=================================


DEALLOCATE(Info%topo%h,Info%topo%hx,Info%topo%hy,Info%topo%hxx,Info%topo%hyy)
DEALLOCATE(Info%topo)
DEALLOCATE(Info%ElliInfo)
DEALLOCATE(Info)
DEALLOCATE(ut,vt,wt)

CALL MKL_FREEBUFFERS()

END

